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statement  of  Problem  Studied 


Algorithms  and  software  for  the  adaptive  solution  of  multiscale  problems  involving  partial 
differential  equations  and  linkage  to  atomic  level  simulations  were  researched.  Techniques  for 
using  the  discontinuous  Galerkin  method  to  solve  hyperbolic  and  singularly  perturbed  parabolic 
problems  were  developed.  New  anisotropic  adaptive  and  parallel  solution  techniques,  a 
posteriori  error  estimation  strategies,  limiting  procedures  that  reduce  spurious  oscillations  near 
discontinuities,  and  discontinuity  detection  strategies  that  reduce  the  need  for  limiting,  thereby 
reducing  both  excess  diffusion  and  spurious  oscillations  were  developed.  The  software  and 
methods  are  being  tested  on  a  variety  of  problems  involving  compressible  flows.  In 
collaboration  with  engineers  at  Benet  Laboratories,  we  have  been  investigating  muzzle  blast 
from  caimons  with  perforated  brakes. 

A  new  procedure  to  couple  atomic/continuum  level  adaptive  simulations  was  developed  and 
demonstrated  on  test  problems.  Scale  error  indicators  have  developed  and  adaptive  construction 
of  local  atomic  regions  demonstrated. 

Summary  of  Most  Important  Results 

Discontinuous  Galerkin  Methods 

An  adaptive  software  system  for  solving  multiscale  hyperbolic  conservation  laws  by  the 
discontinuous  Galerkin  (DG)  method  has  been  developed.  The  software  is  capable  of  addressing 
problems  in  one,  two,  and  three  spatial  dimensions  using  anisotropic  adaptive  /?-refmement. 

The  DG  method  offers  several  advantages  relative  to  traditional  finite  element  and  finite  volume 
methods  when  solving  hyperbolic  systems  of  conservation  laws.  The  advantages  and  several 
aspects  of  the  DG  method  are  described  in  a  sequence  of  manuscripts  that  construct  orthogonal 
bases  [18,21],  describe  limiting  procedures  [10],  develop  efficient  local  time  stepping  algorithms 
[5,18],  create  efficient  serial  and  parallel  data  representations  [20],  and  apply  the  method  to 
compressible  flow  problems  [2,3,18,22]. 

In  order  to  guide  adaptive  enrichment  and  appraise  the  accuracy  of  computed  solutions,  we 
developed  a  posteriori  estimates  of  discretization  errors.  We  show  [1,9]  that  the  leading  term  of 
the  spatial  discretization  error  of  a  piecewise-polynomial  solution  of  degree  p  is  the  difference 
between  orthogonal  polynomials  of  degrees  p  and  p+\.  In  one  dimension,  the  orthogonal 
polynomials  are  Radau  polynomials  and  in  two  dimensions  on  triangular  elements  they  are  linear 
combinations  of  Dubiner  polynomials.  We  further  demonstrate  a  strong  superconvergence  at 
downwind  element  boundaries  where  the  average  spatial  error  converges  as  [9]. 

The  strong  superconvergence  at  outflow  boundaries  can  be  used  to  identify  discontinuities.  This, 
in  turn,  can  be  used  to  provide  information  on  where  to  apply  limiting  to  reduce  oscillations 
when  using  high-order  methods.  This  procedure  has  been  remarkably  successful  at  reducing 
both  excess  dissipation  and  spurious  oscillations  near  discontinuities  [10]. 


To  increase  the  effectiveness  of  the  adaptive  procedures  a  set  of  anisotropic  adaptive  mesh 
control  technologies  have  been  developed.  The  two  key  components  of  this  procedure  are 
anisotropic  correction  indication  procedures  and  anisotropic  mesh  adaptation  procedures.  The 
anisotropic  correction  procedure  has  two  major  components  [22]  that  construct  a  anisotropic 
mesh  metric  tensor  over  the  problem  domain.  In  smooth  portions  the  mesh  metric  tensor  employs 
the  Hessian  of  the  solution  field  [12,13,22].  The  portions  of  the  domain  that  are  not  smooth 
employ  the  shock  capturing  procedure  of  reference  [10].  In  these  locations  the  mesh  metric 
tensor  is  defined  to  align  with  the  shocks  and  to  blend  into  the  smooth  portions  of  the  domain 
[22].  The  anisotropic  mesh  adaptation  is  carried  out  using  a  complete  set  of  mesh  modification 
operations  for  anisotropic  meshes  [12,13]  and  fully  account  for  curved  domain  boundaries  [11]. 


Figures  1  and  2  show  a  2-D  and  3-D  result  obtained  with  the  adaptive  anisotropic  adaptive  DG 
procedure  just  described. 


Figure  1 .  The  meshes  and  pressure  contours  for  t=0.0,  t=0.0002  and  t=0.0005  in  a  2-D  Muzzle  blast 
simulation. 


Figure  2.  Mesh  and  pressure  contours  of  a  3-D  blast  problem. 


Level  Set  Method  for  Solving  Fluid/Rigid  body  Interaction  Problems 

The  simulation  of  flow  problems  with  moving  objects  require  an  effective  means  to  account  for 
the  motion  of  the  object  through  the  domain  of  the  mesh.  A  level  set  based  method  that  allows 
the  object  to  move  through  the  mesh  has  been  developed  to  allow  rigid  objects  to  move  through  a 
mesh.  In  this  method  an  Eulerian  mesh  is  created  over  the  entire  computational  domain  [6].  The 
movement  of  the  rigid  object  will  cut  some  of  the  mesh  entities  and  introduce  new  boundaries 
inside  the  fluid.  The  level  set  technique  is  used  to  track  the  interface  between  the  rigid  body  and 
the  compressible  fluid,  and  the  Ghost  Fluid  Method  (GFM)  to  capture  the  correct  boimdary 
condition  implicitly  at  the  interface.  Therefore,  handling  the  boundaries  is  simplified  and  we 
have  the  advantage  of  doing  a  fluid  calculation  on  the  entire  domain.  By  doing  this,  we  can 
handle  the  contribution  of  the  embedded  moving  boundaries  in  the  fluid  without  mesh 
modification. 

For  the  elements  in  the  Eulerian  mesh,  there  are  two  different  possible  cases  according  to  their 
position  with  respect  to  the  zero  level  set  of  the  level  set  function.  Case  (i),  the  entity  is  totally 
outside  the  zero  level  set.  Case  (ii),  the  entity  is  totally  inside  the  zero  level  set  or  the  entity  is  cut 
by  the  zero  level  set.  For  all  entities  of  case  (i),  the  object  in  the  fluid  has  no  direct  effect  on  this 
entity,  so  no  special  operations  are  needed.  For  all  entities  of  case  (ii),  the  boundary  of  the 
fluid/solid  interface  must  be  accounted  for  using  the  GFM  to  capture  the  correct  boundary 
conditions.  The  level  set  enables  us  to  locate  those  entities  that  need  special  treatment,  even 
when  we  do  not  know  the  position  of  the  precise  position  of  the  moving  boundary  in  the  entity. 
Actually,  we  also  do  not  want  to  find  the  precise  position  of  the  moving  boundary  because  it  is 
an  expensive  operation,  especially  for  high  order  approximation  and  for  high  dimension  case. 

We  treat  the  rigid  body  moving  in  a  fluid  as  a  contact  discontinuity  at  the  interface  because  the 
fluid  near  the  interface  will  move  with  the  velocity  of  the  rigid  body.  The  Rankine-Hugoniot 


jump  conditions  imply  that  both  the  pressure  and  the  normal  veloeity  are  continuous  across  the 
interfaee.  Therefore,  we  set  the  values  for  the  ghost  entity  to  enforce  these  two  eonditions.  There 
are  three  kinds  of  state  variables  that  need  to  be  set:  veloeity,  density,  energy  (or  pressure).  At 
the  same  time,  since  the  interface  of  the  rigid  body  acts  as  a  reflective  wall  boundary  for  the 
fluid,  a  reflective  constraint  can  capture  the  boundary  eondition.  For  DGM,  ghost  values  are  set 
at  every  ghost  entity  by  eonsidering  and  setting  the  values  at  every  gauss  point  inside  the  rigid 
object. 

This  approximation  is  straightforward  even  in  the  multi-dimensional  ease.  This  is  an  advantage 
because  we  ean  treat  the  one-,  two-and  three-dimensional  cases  with  the  same  seheme  and  the 
numerical  implementation  is  simple  and  consistent  for  all  cases.  On  the  other  hand,  a  lot  of 
reflective  points  may  need  to  be  searched.  To  speed  the  procedure,  an  oetree  data  structure  is 
designed  and  implemented.  In  the  case  that  a  reflective  point  does  not  belong  to  the 
computational  domain  of  the  problem,  specific  care  must  be  taken  in  setting  values. 

This  procedure  has  been  eombined  with  adaptive  mesh  refinement  [21,22]  to  eontrol  solution 
errors  to  solve  problems  with  moving  objeets  (see  Figure  3). 


Figure  3.  Example  sphere  moving  through  a  fluid  domain. 

Coupled  Atomistic-Continuum  Simulations 

We  have  developed  a  Composite  Atomistic-Continuum  Method  (CACM)  for  multiscale  analysis 
of  eoupled  discrete  and  continuum  models  [5],  The  main  idea  of  the  method  is  to  use  atomistic 
analysis  only  at  places  where  needed  to  capture  the  highly  nonlinear  behavior  and  continuum 
analysis  will  be  used  elsewhere  for  efficieney.  Such  a  method  permits  the  explicit  representation 
of  the  nano-scale  physies.  The  components  of  the  CACM  and  its  development  include: 

1.  CACM  Method:  In  CACM,  the  continuum  is  represented  by  linear  elastieity  model,  and 
diseretized  using  finite  elements  on  a  coarse  grid.  An  atomistie  representation  in  the 
regions  of  high  field  gradients  or  where  the  highly  nonlinear  and  non-loeal  behavior  of 
the  lattice  is  important,  i.e.,  at  the  places  where  there  is  nucleation  or  close  range 
interaction  of  the  defects  is  used.  In  these  regions,  an  atomistic  grid  is  used  which  carries 
the  boundary  conditions  imposed  by  the  continuum  and  returns  to  the  continuum 
correcting  information.  The  atomistic  grid  is  defined  in  restricted  spatial  regions  and 


relocates  with  the  defects.  This  procedure  minimizes  the  computational  effort  associated 
with  atomistic  simulations,  while  insuring  the  required  accuracy  at  relevant  spatial 
locations.  We  have  tried  various  algorithms,  and  operators  to  link  the  two  scales.  Based 
on  the  effectiveness,  meaningfulness  and  ease  of  transferring  the  information  from  the 
continuum  to  the  atomistic,  and  back  to  the  continuum,  at  this  point  we  have  selected  to 
transfer  the  displacement  information  from  the  continuum  to  the  atomistic  as  the 
boundary  condition  for  the  atomistic  analysis.  Then  in  turn  the  atomistic  analysis 
provides  the  correction  of  the  continuum  solution  at  the  critical  places.  The  atomistic  grid 
is  overlaid  on  the  continuum  grid,  and  this  method  doesn't  have  the  restriction  of 
matching  the  two  grids.  This  allows  for  using  quite  coarse  grid  at  the  continuum  level 
unlike  some  other  methods,  which  is  forced  to  use  refined  meshes  for  the  continuum  near 
the  atomistic  regions  to  match  the  atomistic  grid. 

2.  Implementation  of  the  CACM:  The  algorithm  for  the  CACM  was  implemented  in  the 
environment  of  Trellis,  a  framework  for  3D  adaptive  multi -physics  and  multi-scale 
analysis.  The  Trellis  framework  was  extended  to  build  the  scope  for  multiscale 
applications.  The  capability  to  solve  3D  atomistic  problems  using  various  types  of 
potentials  was  implemented  within  the  Trellis  framework.  To  couple  the  continuum 
analysis  with  the  atomistic  analysis  a  module  responsible  for  transferring  the  information 
between  the  two  analyses  was  implemented.  This  module  has  been  implemented  in  a 
general  manner  so  as  to  be  able  to  be  utilized  for  possible  coupling  of  other  analyses. 

3.  Error  indicators:  We  have  identified  the  following  three  criterions  required  for  the 
adaptive  selection  of  the  model: 

(a)  Non  linearity:  This  is  characterized  by  the  magnitude  of  the  strain.  We  have  already 
developed  an  error  indicator  to  capture  this  property.  The  error  indicator  is  based  on 
the  relative  magnitude  of  the  first  order  energy  term  in  comparison  to  the  zeroth  order 
term. 

(b)  Non  locality:  This  is  characterized  by  the  gradient  of  the  strain.  Indicators  to  capture 
this  information  were  proposed,  and  used  effectively  in  their  method. 

(c)  Non-convexity:  This  is  characterized  by  the  change  in  energy  felt  by  the  continuum 
once  a  defect  passes  through  it. 

4.  Adaptive  modeling:  The  software  has  to  be  generalized  in  the  following  respeets  to 
handle  adaptivity:  (a)  Building  a  module  that  will  provide  indicators  of  the  error  in 
modeling,  (b)  Automatic  generation  of  the  discrete  regions  based  on  the  error 
information. 

5.  Demonstration  of  the  CACM:  Figure  4  shows  an  example  of  the  multiscale  modeling  of 
the  growth  of  defects  at  a  crack  tip  [5], 

Technology  Transfer 

Joseph  E.  Flaherty  and  Mark  S.  Shephard  have  regular  interactions  with  Robert  Dillon,  Deborah 
Bleau,  and  Daniel  Cler  of  Benet  Laboratories  [2,3].  As  described  in  Section  1,  we  have  been 
eollaborating  on  using  the  DG  software  to  analyze  muzzle  blast  effects  with  large-calibre 
weapons  systems.  These  problems  involve  complex  three-dimensional  interior  and  exterior 
flows  including  realistic  geometries  and  ground  reflections. 


a)3-D  domain  with  a  crack  and  macroscale  mesh  b)  Adaptively  superimposed  atomic  region  (in  green) 


Figure  4.  Example  of  the  multiscale  modeling  of  the  growth  of  defects  at  a  crack  tip. 

J.  E.  Flaherty  was  also  the  dissertation  advisor  (with  Henry  Nagamatsu)  of  Eric  Kathe,  a 
mechanical  engineer  at  Benet  Laboratories.  Dr.  Kathe  completed  his  Ph.D.  dissertation  in  May 
2002. 

.  A  start-up  company,  Simmetrix  Corporation,  is  building  on  the  mesh  modification  and  analysis 
framework  developments  done  in  this  project.  Simmetrix  has  commercialized  a  number  of  the 
mesh  generation  and  mesh  adaptation  procedures  developed  in  this  project.  They  have  also 
obtained  two  ARMY  supported  SBIR’s.  The  first  is  concerned  with  the  development  of 
commercial  parallel  adaptive  simulation  technologies.  The  second,  which  has  just  recently 
started  is  concerned  with  the  development  of  multiscale  simulation  technologies  with  an 


emphasis  on  atomic  level  simulation  model  generation  and  control.  Rensselaer  is  a  sub¬ 
contractor  to  Simmetrix  on  both  of  these  SBIR’s. 
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